function [h0Adj,hxAdj,sigmaAdj] = bootstrapVAR(factors,h0OLS,hxOLS,sigmaOLS,numBoot,seedNum)

% Resampling
rng(seedNum,'twister')
[nx,T]   = size(factors);

% The residuals
residuals       = zeros(nx,T);
for t=2:T
    residuals(:,t) =  factors(:,t) - (h0OLS + hxOLS*factors(:,t-1));
end

% Setting dimensions
h0Boot    = zeros(nx,numBoot);
hxBoot    = zeros(nx,nx,numBoot);
sigmaBoot = zeros(nx,nx,numBoot);
for s=1:numBoot
    indexBoot         = randi(T,1,T);
    factors_mc        = zeros(nx,T);
    factors_mc(:,1)   = factors(:,indexBoot(1,1));
    for t=2:T
        factors_mc(:,t)   = h0OLS + hxOLS*factors_mc(:,t-1) + residuals(:,indexBoot(1,t));
    end
    
    [B,OMEGA] = Var1(factors_mc');
    h0Boot(:,s)      = B(:,nx+1);
    hxBoot(:,:,s)    = B(1:nx,1:nx);
    sigmaBoot(:,:,s) = chol(OMEGA,'lower');
end

% Average in bootstrap
meanh0Boot = mean(h0Boot,2);
meanhxBoot = mean(hxBoot,3);
meansigmaBoot = mean(sigmaBoot,3);

% The bias adjusted estimates
h0Adj     = h0OLS    - (meanh0Boot - h0OLS);
hxAdj     = hxOLS    - (meanhxBoot - hxOLS);
sigmaAdj  = sigmaOLS - (meansigmaBoot - sigmaOLS);



end

